# ------------------------------------------------------------------------------
##### Set up                                                               -----
# ------------------------------------------------------------------------------

##### Libraries ----------------------------------------------------------------
library(tidyverse)
library(stargazer)
library(marginaleffects)
library(readr)
library(ggtext)
library(ggrepel)
library(broom)

# Function to get summary stats
get_descriptives <- function(data, variables) {
  descriptives <- lapply(variables, function(var) {
    x <- data[[var]]
    c(
      min = min(x, na.rm = TRUE),
      max = max(x, na.rm = TRUE),
      median = median(x, na.rm = TRUE),
      mean = mean(x, na.rm = TRUE),
      sd = sd(x, na.rm = TRUE),
      n = sum(!is.na(x)),
      missing = sum(is.na(x))
    )
  })
  
  result <- do.call(rbind, descriptives)
  rownames(result) <- variables
  as.data.frame(result)
}

##### Load data ----------------------------------------------------------------
difdef <- read_csv("data/difdef.csv")

# ------------------------------------------------------------------------------
##### Descriptives                                                         -----
# ------------------------------------------------------------------------------

# N per treatmentgroup                                                     -----
n_actor <- difdef |> 
  group_by(actor) |> 
  summarise(count = n())

stargazer(n_actor,
          summary = FALSE,
          type = "text",
          out = "tables/n_actor.tex",
          label = "tab:n_actor",
          title = "N per specific actor")

n_type <- difdef |> 
  group_by(type) |> 
  summarise(count = n())

stargazer(n_type,
          summary = FALSE,
          type = "text",
          out = "tables/n_type.tex",
          label = "tab:n_type",
          title = "N per type of democratic defender")

# DVs                                                                      -----
dvs <- c("eval_pre", "eval_post", "outcome",
         "perception_neut", "perception_recog", "perception_leri", "perception_comt",
         "leri_bias")

summary_dvs <- get_descriptives(difdef, dvs)

stargazer(summary_dvs,
          summary = FALSE,
          type = "text",
          out = "tables/desc_dvs.tex",
          label = "tab:desc_dvs",
          title = "Descriptives for the outcome variables")

# IVs (interval)                                                           -----
ivs <- c("leri", "dem_satis", "dem_sup", "pol_interest", "pol_news", "pol_trust_inparty",
         "pol_trust_outparty", "pol_trust_judges", "pol_trust_govt", "pol_trust_parl",
         "pol_trust_acad", "pol_trust_ngo", "pol_trust_mili")

summary_ivs <- get_descriptives(difdef, ivs)

stargazer(summary_ivs,
          summary = FALSE,
          type = "text",
          out = "tables/desc_ivs.tex",
          label = "tab:desc_ivs",
          title = "Descriptives for the independent variables (interval)")

# IVs (categorical)                                                        -----
employment <- difdef |> 
  mutate(employment = case_when(employment %in% c("CONSENT_REVOKED", "DATA_EXPIRED") ~ "No data", # This is respondents' consent to Prolific to share demographics, not respondents' consent to participate in the survey.
                                TRUE ~ employment)) |> 
  group_by(employment) |> 
  summarise(count = n())

stargazer(employment,
          summary = FALSE,
          type = "text",
          out = "tables/desc_employment.tex",
          label = "tab:desc_employment",
          title = "Descriptives: Employment")

# Sample representativeness                                                -----
# population_counts provided by Prolific,
# more information https://researcher-help.prolific.com/en/article/e6555f#:~:text=Where%20does%20Prolific%20get%20its%20census%20data%20from%3F
# archived pdf of the webpage is available in the dataverse

# UK census
# England and Wales (2021): https://www.nomisweb.co.uk/datasets/c2021rm032
# Northern Ireland (2021): https://www.nisra.gov.uk/publications/census-2021-main-statistics-ethnicity-tables
# Scotland (2011): https://www.scotlandscensus.gov.uk/documents/2011-census-table-data-all-areas/
age <- difdef |> 
  group_by(Age) |> 
  summarise(sample_count = n()) |> 
  mutate(sample_percentage = round(sample_count / sum(sample_count) * 100, 1))

age_repr <- tibble(
  Age = c("18-24", "25-34", "35-44", "45-54", "55+"),
  pop_count = c(422, 671, 648, 664, 1567)) |> 
  mutate(pop_percentage = round(pop_count / sum(pop_count) * 100, 1)) |> 
  left_join(age, by = "Age") |> 
  dplyr::select(Age, sample_count, sample_percentage, pop_percentage) |> 
  rename("Sample count" = sample_count,
         "Sample %" = sample_percentage,
         "Population %" = pop_percentage)

stargazer(age_repr,
          summary = FALSE,
          type = "text",
          out = "tables/desc_age.tex",
          label = "tab:desc_age",
          title = "Descriptives: Age")

sex <- difdef |> 
  mutate(Sex = case_when(Sex == "CONSENT_REVOKED" ~ "No data", # This is respondents' consent to Prolific to share demographics, not respondents' consent to participate in the survey.
                         TRUE ~ Sex)) |>  
  group_by(Sex) |> 
  summarise(sample_count = n()) |> 
  mutate(sample_percentage = round(sample_count / sum(sample_count) * 100, 1))

sex_repr <- tibble(
  Sex = c("Female", "Male", "No data"),
  pop_count = c(2051, 1921, 0)) |> 
  mutate(pop_percentage = round(pop_count / sum(pop_count) * 100, 1)) |> 
  left_join(sex, by = "Sex") |> 
  dplyr::select(Sex, sample_count, sample_percentage, pop_percentage) |> 
  rename("Sample count" = sample_count,
         "Sample %" = sample_percentage,
         "Population %" = pop_percentage)

stargazer(sex_repr,
          summary = FALSE,
          type = "text",
          out = "tables/desc_sex.tex",
          label = "tab:desc_sex",
          title = "Descriptives: Sex")

ethnicity <- difdef |> 
  mutate(ethnicity = case_when(ethnicity == "CONSENT_REVOKED" ~ "No data", # This is respondents' consent to Prolific to share demographics, not respondents' consent to participate in the survey.
                         TRUE ~ ethnicity)) |> 
  group_by(ethnicity) |> 
  summarise(sample_count = n()) |> 
  mutate(sample_percentage = round(sample_count / sum(sample_count) * 100, 1))

etn_repr <- tibble(
  ethnicity = c("Asian", "Black", "Mixed", "No data", "Other", "White"),
  pop_count = c(312, 131, 66, 0, 72, 3391)) |> 
  mutate(pop_percentage = round(pop_count / sum(pop_count) * 100, 1)) |> 
  left_join(ethnicity, by = "ethnicity") |> 
  dplyr::select(ethnicity, sample_count, sample_percentage, pop_percentage) |> 
  rename("Ethnicity" = ethnicity,
         "Sample count" = sample_count,
         "Sample %" = sample_percentage,
         "Population %" = pop_percentage)

stargazer(etn_repr,
          summary = FALSE,
          type = "text",
          out = "tables/desc_ethnicity.tex",
          label = "tab:desc_ethnicity",
          title = "Descriptives: Ethnicity")

##### Expected democratic defenders --------------------------------------------
expected <- read_csv("data/expected.csv")

expected_count <- expected |>
  group_by(ResponseId) |>
  summarise(count = n())

avg_expectation <- mean(expected_count$count)
sd_expectation <- sd(expected_count$count)

no_expectation <- difdef |>
  filter(expected == "no")

actor_counts <- difdef |>
  group_by(actor) |>
  summarise(actor_count = n())

expected_count <- expected |>
  group_by(expected) |>
  summarise(count = n()) |>
  mutate(perc = count / nrow(difdef) * 100,
         label = paste0("*", round(perc, 2), "%*")) |>
  filter(!is.na(expected)) |>
  mutate(name = "expected") |>
  arrange(count) |>
  dplyr::select(-name)

# plot parameters
expected_plot <- bind_rows(expected_count) |>
  mutate(label_y = perc + 1.5,
         expected = case_when(expected == "judge" ~ "Judge",
                              expected == "civil servant" ~ "Civil servant",
                              expected == "military official" ~ "Military official",
                              expected == "out-party" ~ "Out-party politician",
                              expected == "in-party" ~ "In-party politician",
                              expected == "mayor" ~ "Local politician",
                              expected == "academic" ~ "Academic",
                              expected == "journalist" ~ "Journalist",
                              expected == "NGO" ~ "International  \norganization",
                              expected == "other" ~ "Other",
                              TRUE ~ expected))


# set order
expected_plot$expected <- factor(expected_plot$expected, levels = c("International  \norganization", "Journalist", "Academic",
                                                                    "Local politician", "In-party politician", "Out-party politician",
                                                                    "Military official", "Civil servant", "Judge",
                                                                    "Other"))

expected_colour <- c("expected" = "#C33C54",
                     "inline"= "#37718E")

expected_shape <- c("expected" = 22,
                    "inline"= 23)

expected_labels <- c("expected" = "Percentage of respondents who  \nexpect the actor to defend democracy",
                     "inline" = "Percentage of respondents who  \nsaw the actor in the experiment  \nand expect them to defend democracy")

# plot
ggplot(data = expected_plot,
       aes(x = reorder(expected, perc),
           y = perc)) +

  geom_point(colour = "#C33C54",
             fill = "#C33C54",
             shape = 22,
             size = 2,
             position = position_dodge(width = 0.5)) +

  geom_segment(aes(x = expected,
                   y = 0,
                   yend = perc),
               linewidth = 0.5,
               colour = "#C33C54",
               position = position_dodge(width = 0.5)) +
  geom_richtext(aes(label = label,
                    y = label_y),
                colour = "#C33C54",
                fill = NA,
                label.colour = NA,
                size = 2,
                #nudge_y = 1.5,
                position = position_dodge(width = 0.5)) +
  theme_bw() +
  labs(x = "",
       y = "") +
  theme(axis.text.x = element_markdown(angle = 90, hjust = 1),
        legend.position = "none",
        legend.title = element_blank(),
        legend.text = element_markdown(),
        legend.text.position = "bottom",
        legend.key.spacing.x = unit(1.5, "cm"),
        legend.location = "plot",
        axis.ticks.x = element_blank(),
        axis.line = element_line(colour = "#55505C"),
        axis.ticks.y = element_line(colour = "#55505C")) +
  scale_y_continuous(limits = c(0, 55), expand = c(0, 0))

# and save
ggsave(filename  = "figures/expected_democratic_defenders.png",
       width = 16,
       height = 12,
       dpi = 300,
       units = "cm")

expected_table <- expected_count |>
  dplyr::select(-label) |>
  mutate(perc = round(perc, 2)) |>
  rename(Actor = expected,
         Count = count,
         "% of respondents" = perc) |>
  as_tibble()

stargazer(expected_table,
          summary = FALSE,
          type = "latex",
          out = "tables/expected.tex",
          label = "tab:expected",
          title = "Which elites are expected to defend democracy?")

# "Other" expected actor hand-coded -----
other <- read_excel("data/other-democratic-defenders-coded.xlsx") |>
  mutate(code = str_sub(Code, 1, 1),
         defender = case_when(code == "1" ~ "No one",
                              code == "2" ~ "Citizens",
                              code == "3" ~ "Everyone",
                              code == "4" ~ "People",
                              code == "5" ~ "Voters",
                              code == "6" ~ "Specific actors",
                              code == "7" ~ "Other comments",
                              code == "8" ~ "Don't know",
                              code == "9" ~ "Uninformative input"),
         extra = case_when(Code %in% c("1a", "6g", "7d", "8a") ~ "Disillusionment",
                           Code %in% c("3b", "4a") ~ "People vs elite"))

other2 <- other |>
  group_by(defender) |>
  summarise(count = n()) |>
  arrange(desc(count))

stargazer(other2,
          summary = FALSE,
          type = "latex",
          out = "tables/other-democratic-defenders.tex",
          label = "tab:other1",
          title = "Other expected democratic defenders")

other3 <- other |>
  group_by(extra) |>
  filter(!is.na(extra)) |>
  summarise(count = n()) |>
  arrange(desc(count))

stargazer(other3,
          summary = FALSE,
          type = "latex",
          out = "tables/other-extra.tex",
          label = "tab:other2",
          title = "Expressed attitude towards democracy")

# ------------------------------------------------------------------------------
##### Deviation from the preregistration: different leri-measure           -----
# ------------------------------------------------------------------------------
# Results not interpreted in main text
# See Online Appendix

m_dist <- lm(leri_distance ~ actor + inparty + mayor + opposite,
             data = difdef)

dist_pred <- avg_predictions(m_dist, by = "actor") |> 
  mutate(perception = "prereg",
         iv = "actor")

m_dist_type <- lm(leri_distance ~ type + inparty + mayor + opposite,
                  data = difdef)

dist_type_pred <- avg_predictions(m_dist_type, by = "type") |> 
  mutate(perception = "prereg",
         iv = "type")

leri_prereg <- bind_rows(dist_pred, dist_type_pred, leri_pred, leri_type_pred) |> 
  mutate(type = case_when(type == "institutional" ~ "Institutional",
                          type == "social" ~ "Social",
                          type == "political" ~ "Political"),
         actor = ifelse(is.na(actor), type, actor),
         actor = case_when(actor == "institutional" ~ "Institutional",
                           actor == "political" ~ "Political",
                           actor == "social" ~ "Social",
                           actor == "judge" ~ "Judge",
                           actor == "civil servant" ~ "Civil servant",
                           actor == "military official" ~ "Military official",
                           actor == "out-party" ~ "Out-party politician",
                           actor == "in-party" ~ "In-party politician",
                           actor == "mayor" ~ "Mayor",
                           actor == "academic" ~ "Academic",
                           actor == "journalist" ~ "Journalist",
                           actor == "NGO" ~ "International human  \nrights watchdog",
                           TRUE ~ actor),
         label = paste0(round(estimate, 2), " [", round(conf.low, 2), "; ", round(conf.high, 2), "]"))

### plot parameters
# set order
leri_prereg$actor <- factor(leri_prereg$actor, levels = c("International human  \nrights watchdog", "Journalist", "Academic",
                                                          "Mayor", "In-party politician", "Out-party politician",
                                                          "Military official", "Civil servant", "Judge",
                                                          "Social", "Political", "Institutional"))

leri_prereg$perception <- factor(leri_prereg$perception, levels = c("left-right", "prereg"))

leri_prereg$iv <- factor(leri_prereg$iv, levels = c("type", "actor"))

# colours
perception_colours <- c("prereg" = "#C33C54",
                        "left-right" = "#37123C")

perception_labels <- c("prereg" = "**Left-right bias as distance**  \n(preregistered, but omitted in main text)",
                       "left-right" = "**Folded left-right bias**  \n(interpreted in main text)") 

perception_shapes <- c("prereg" = 22,
                       "left-right" = 23)

# plot
ggplot(data = leri_prereg,
       aes(x = estimate,
           y = actor)) +
  
  geom_errorbarh(aes(xmin = conf.low, 
                     xmax = conf.high,
                     colour = perception),
                 height = 0,
                 position = position_dodge(width = 0.5),
                 linewidth = 0.2) +
  
  geom_point(aes(fill = perception,
                 shape = perception,
                 colour = perception),
             position = position_dodge(width = 0.5),
             size = 1) +
  
  geom_text(aes(label = label,
                colour = perception),
            size = 2.5,
            position = position_dodge(width = 1.3),
            show.legend = FALSE) +
  
  facet_grid(iv ~ .,
             scales = "free_y",
             space = "free",
             labeller = labeller(iv = c("actor" = "Specific actor", "type" = "Category"))) +
  
  # formatting
  theme_bw() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_markdown(),
        legend.text.position = "bottom",
        legend.key.spacing.x = unit(1.5, "cm"),
        legend.location = "plot",
        panel.grid.major.x = element_line(),
        axis.text.y = element_markdown(),
        axis.title = element_blank(),
        strip.text = element_markdown(),
        axis.ticks.y = element_blank()) +
  scale_x_continuous(limits = c(1, 5), breaks = 1:8) +
  scale_fill_manual(values = perception_colours, labels = perception_labels) +
  scale_colour_manual(values = perception_colours, labels = perception_labels) +
  scale_shape_manual(values = perception_shapes, labels = perception_labels)

# save
ggsave(filename  = "figures/perceptions_prereg.png",
       width = 20,
       height = 16,
       dpi = 300,
       units = "cm")

# ------------------------------------------------------------------------------
##### Deviation from the preregistration: results for omitted H2           -----
# ------------------------------------------------------------------------------
# Results not interpreted in main text
# See Online Appendix

##### Difference (within-subject) ----------------------------------------------
h2a <- lm(outcome ~ perception_neut + inparty + mayor + opposite,
          data = difdef)

h2b <- lm(outcome ~ leri_bias + inparty + mayor + opposite,
          data = difdef)

h2c <- lm(outcome ~ perception_recog + inparty + mayor + opposite,
          data = difdef)

h2d <- lm(outcome ~ perception_comt + inparty + mayor + opposite,
          data = difdef)


# Table
stargazer(h2c, h2d,
          title = "Effect of embedded knowledge-perceptions on democracy evaluation (within individuals)",
          column.labels = c("Ability to recognize", "Commitment to democracy"),
          dep.var.caption = "",
          dep.var.labels.include = FALSE,
          align = TRUE,
          omit.stat = c("LL", "ser", "f"),
          no.space = TRUE,
          #single.row = TRUE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type = "latex",
          out = "tables/h2_ek.tex",
          label = "tab:h2_ek")

stargazer(h2a, h2b,
          title = "Effect of political neutrality-perceptions on democracy evaluation (within individuals)",
          column.labels = c("Partisanship", "Left-right distance"),
          dep.var.caption = "",
          dep.var.labels.include = FALSE,
          align = TRUE,
          omit.stat = c("LL", "ser", "f"),
          no.space = TRUE,
          #single.row = TRUE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type = "latex",
          out = "tables/h2_pn.tex",
          label = "tab:h2_pn")

##### Post-only (between-group) ------------------------------------------------

h2a_post <- lm(eval_post ~ perception_neut + inparty,
               data = difdef)

h2b_post <- lm(eval_post ~ leri_bias + inparty,
               data = difdef)

h2c_post <- lm(eval_post ~ perception_recog + inparty,
               data = difdef)

h2d_post <- lm(eval_post ~ perception_comt + inparty,
               data = difdef)

stargazer(h2c_post, h2d_post,
          title = "Effect of embedded knowledge-perceptions on democracy evaluation (between groups)",
          column.labels = c("Ability to recognize",  "Commitment to democracy"),
          dep.var.caption = "",
          dep.var.labels.include = FALSE,
          align = TRUE,
          omit.stat = c("LL", "ser", "f"),
          no.space = TRUE,
          #single.row = TRUE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type = "latex",
          out = "tables/h2_ek_post.tex",
          label = "tab:h2_ek_post")

stargazer(h2a_post, h2b_post,
          title = "Effect of political neutrality-perceptions on democracy evaluation (between groups)",
          column.labels = c("Partisanship", "Left-right distance"),
          dep.var.caption = "",
          dep.var.labels.include = FALSE,
          align = TRUE,
          omit.stat = c("LL", "ser", "f"),
          no.space = TRUE,
          #single.row = TRUE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type = "latex",
          out = "tables/h2_pn_post.tex",
          label = "tab:h2_pn_post")

# Prepare plot data
h2a_plot <- summary(h2a) |> 
  tidy() |> 
  filter(term == "perception_neut") |> 
  mutate(model = "diff")

h2b_plot <- summary(h2b) |> 
  tidy() |> 
  filter(term == "leri_bias") |> 
  mutate(model = "diff")

h2c_plot <- summary(h2c) |> 
  tidy() |> 
  filter(term == "perception_recog")|> 
  mutate(model = "diff")

h2d_plot <- summary(h2d) |> 
  tidy() |> 
  filter(term == "perception_comt")|> 
  mutate(model = "diff")

h2a_post_plot <- summary(h2a_post) |> 
  tidy() |> 
  filter(term == "perception_neut") |> 
  mutate(model = "post")

h2b_post_plot <- summary(h2b_post) |> 
  tidy() |> 
  filter(term == "leri_bias") |> 
  mutate(model = "post")

h2c_post_plot <- summary(h2c_post) |> 
  tidy() |> 
  filter(term == "perception_recog")|> 
  mutate(model = "post")

h2d_post_plot <- summary(h2d_post) |> 
  tidy() |> 
  filter(term == "perception_comt")|> 
  mutate(model = "post")

h2_plot <- bind_rows(h2a_plot, h2b_plot, h2c_plot, h2d_plot, 
                     h2a_post_plot, h2b_post_plot, h2c_post_plot, h2d_post_plot) |> 
  mutate(order = case_when(term == "perception_recog" ~ 1,
                           term == "perception_comt" ~ 2,
                           term == "perception_neut" ~ 3,
                           term == "leri_bias" ~ 4),
         term = case_when(term == "perception_recog" ~ "Perceived ability to\n recognize autocratization",
                          term == "perception_comt" ~ "Perceived commitment\nto democracy",
                          term == "perception_neut" ~ "Perceived partiality",
                          term == "leri_bias" ~ "Perceived left-right bias"),
         label = paste0(round(estimate, 2), " [", round(estimate - 1.96 * std.error, 2), "; ", round(estimate + 1.96 * std.error, 2), "]"))


# Plot parameters
h2_colour <- c("diff" = "#C33C54",
               "post"= "#37718E")

h2_shape <- c("diff" = 22,
              "post"= 23)

h2_labels <- c("diff" = "**Primary outcome**  \n_(Within-subject difference between  \npre- and post-treatment evaluation)_",
               "post" = "**Secondary outcome**  \n_(Between-group post-treatment evaluation)_") 

# Plot
ggplot(data = h2_plot,
       aes(x = estimate,
           y = reorder(term, order))) +
  geom_vline(xintercept = 0,
             colour = "darkgrey",
             linewidth = 0.2) +
  geom_point(aes(colour = model,
                 shape = model,
                 fill = model),
             position = position_dodge2(width = 0.5)) +
  geom_errorbarh(aes(xmin = estimate - 1.96 * std.error,
                     xmax = estimate + 1.96 * std.error,
                     colour = model),
                 position = position_dodge(width = 0.5),
                 height = 0) +
  geom_text(aes(label = label),
            position = position_dodge2(width = 1),
            size = 2) +
  theme_bw() +
  labs(x = "",
       y = "") +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_markdown(),
        legend.text.position = "bottom",
        legend.key.spacing.x = unit(1.5, "cm"),
        legend.location = "plot",
        axis.title.x = element_markdown(),
        axis.ticks.x = element_blank(),
        axis.line = element_line(colour = "#55505C"),
        axis.ticks.y = element_line(colour = "#55505C")) +
  scale_colour_manual(values = h2_colour, labels = h2_labels) +
  scale_fill_manual(values = h2_colour, labels = h2_labels) +
  scale_shape_manual(values = h2_shape, labels = h2_labels)

# ... and save!
ggsave(filename  = "figures/h2_diff_and_post.png",
       width = 18,
       height = 12,
       dpi = 300,
       units = "cm")

# # ------------------------------------------------------------------------------
# ##### Hypothesis 1                                                         -----
# # ------------------------------------------------------------------------------
# h1_sub <- lm(outcome ~ type + inparty + mayor + opposite,
#          data = difdef |> filter(eval_pre > 1))
# 
# h1_sub_pred <- avg_predictions(h1_sub, by = "type") |> 
#   rename(actor = type) |> 
#   mutate(model = "agg")
# 
# h1_actor_sub <- lm(outcome ~ actor + inparty + mayor + opposite,
#                data = difdef |> filter(eval_pre > 1))
# 
# h1_actor_sub_pred <- avg_predictions(h1_actor_sub, by = "actor") |> 
#   mutate(model = "dis")
# 
# h1_sub_plot <- bind_rows(h1_sub_pred, h1_actor_sub_pred) |> 
#   mutate(stars = case_when(p.value < 0.05 & p.value >= 0.01 ~ "*",
#                            p.value < 0.01 & p.value >= 0.001 ~ "**",
#                            p.value < 0.001 ~ "***",
#                            TRUE ~ ""),
#          label = paste0(round(estimate, 2), stars, " [", round(conf.low, 2), "; ", round(conf.high, 2), "]")) |> 
#   mutate(actor = case_when(actor == "institutional" ~ "**Institutional**",
#                            actor == "political" ~ "**Political**",
#                            actor == "social" ~ "**Social**",
#                            TRUE ~ actor))
# 
# # set order
# h1_sub_plot$actor <- factor(h1_sub_plot$actor, levels = c("NGO", "journalist", "academic", "**Social**",
#                                                   "mayor", "in-party", "out-party", "**Political**",
#                                                   "military official", "civil servant", "judge", "**Institutional**"))
# 
# # plot parameters
# h1_colours <- c("agg" = "#C33C54",
#                 "dis" = "#37718E")
# 
# h1_fill <- c("agg" = "#C33C54",
#              "dis" = "white")
# 
# ggplot(data = h1_sub_plot,
#        aes(x = estimate,
#            y = actor)) +
#   geom_vline(xintercept = 0,
#              colour = "darkgrey",
#              linewidth = 0.2) +
#   geom_errorbarh(aes(xmin = conf.low,
#                      xmax = conf.high,
#                      colour = model),
#                  height = 0) +
#   geom_point(aes(colour = model,
#                  fill= model),
#              size = 3,
#              shape = 21) +
#   geom_text(aes(label = label),
#             nudge_y = 0.3,
#             size = 3) +
#   scale_colour_manual(values = h1_colours) +
#   scale_fill_manual(values = h1_fill) +
#   theme_classic() +
#   theme(legend.position = "none",
#         axis.ticks.y = element_blank(),
#         axis.text.y = element_markdown()) +
#   labs(x = "Effectiveness of democratic defence",
#        y = "")
# 
# ggsave(filename  = "figures/h1.png",
#        width = 16,
#        height = 12,
#        dpi = 300,
#        units = "cm")
# 
# stargazer(h1, h1_actor,
#           title = "Effect of democratic defenders on democracy evaluation",
#           column.labels = c("Aggregate type of actor", "Specific actor"),
#           align = TRUE,
#           omit.stat = c("LL", "ser", "f"),
#           no.space = TRUE,
#           single.row = TRUE,
#           star.cutoffs = c(0.05, 0.01, 0.001),
#           type = "text",
#           out = "tables/h1.html")

# ------------------------------------------------------------------------------
##### Hypothesis 1  (post-only, between-group)                             -----
# ------------------------------------------------------------------------------

h1_post <- lm(eval_post ~ type + inparty + mayor + opposite,
              data = difdef)

h1_post_actor <- lm(eval_post ~ actor + inparty + mayor + opposite,
                    data = difdef)

stargazer(h1_post, h1_post_actor,
          title = "Effect of democratic defenders on democracy evaluation",
          column.labels = c("Aggregate type of actor", "Specific actor"),
          dep.var.caption = "",
          dep.var.labels.include = FALSE,
          align = TRUE,
          omit.stat = c("LL", "ser", "f"),
          no.space = TRUE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          keep = c("type.*", "actor"),
          notes = "Some variables omitted from table for clarity",
          type = "latex",
          out = "tables/h1_post.tex",
          label = "tab:h1_post")

# ------------------------------------------------------------------------------
##### Attention and Manipulation                                           -----
# ------------------------------------------------------------------------------
attman <- difdef |> 
  dplyr::select(ResponseId, actor, type, attention, manipulation,
                outcome, starts_with("perception_"))

# Attention and manipulation counts
att <- attman |> 
  group_by(attention) |> 
  summarise(count = n()) |> 
  mutate(percentage = paste0(round(100 * count / sum(count), 2), "%"))

stargazer(att,
          type = "latex",
          out = "tables/attention.tex",
          title = "Attention",
          label = "tab:att",
          summary = FALSE)

man <- attman |> 
  group_by(manipulation) |> 
  summarise(count = n()) |> 
  mutate(percentage = paste0(round(100 * count / sum(count), 2), "%"))

stargazer(man,
          type = "latex",
          out = "tables/manipulation.tex",
          title = "Manipulation",
          label = "tab:man",
          summary = FALSE)

# ------------------------------------------------------------------------------
##### Heterogeneous Treatment Effects                                      -----
# ------------------------------------------------------------------------------

# Define models for HTEs
models <- list(
  h1 = outcome ~ type + inparty + mayor + opposite,
  h2_neutrality = perception_neut ~ type + inparty + mayor + opposite,
  h2_bias = leri_bias ~ type + inparty + mayor + opposite,
  h3_recognition = perception_recog ~ type + inparty + mayor + opposite,
  h3_commitment = perception_comt ~ type + inparty + mayor + opposite)

# Different models excluding in-party for party-HTEs
party_models <- list(
  h1 = outcome ~ type + mayor + opposite,
  h2_neutrality = perception_neut ~ type + mayor + opposite,
  h2_bias = leri_bias ~ type + mayor + opposite,
  h3_recognition = perception_recog ~ type + mayor + opposite,
  h3_commitment = perception_comt ~ type + mayor + opposite)

# LeRi -------------------------------------------------------------------------
leri_subgroups <- list(
  left = filter(difdef, leri <= 3),
  center = filter(difdef, leri >= 4 & leri <= 6),
  right = filter(difdef, leri >= 7)
)

leri_labels <- c("Left-leaning respondents", "Center respondents", "Right-leaning respondents")

run_leri_models <- function(model_name, formula) {
  
  fits <- lapply(leri_subgroups, function(data) lm(formula, data = data))
  
  stargazer(fits[[1]], fits[[2]], fits[[3]],
            title = paste0(toupper(model_name), " for left-right sub-groups"),
            column.labels = leri_labels,
            dep.var.caption = "",
            dep.var.labels.include = FALSE,
            align = TRUE,
            omit.stat = c("LL", "ser", "f"),
            no.space = TRUE,
            star.cutoffs = c(0.05, 0.01, 0.001),
            keep = c("type.*", "perception.*", "leri.*", "(Intercept)"),
            notes = "Some variables omitted from table for clarity",
            type = "latex",
            out = paste0("tables/", model_name, "_hte_leri.tex"),
            label = paste0("tab:", model_name, "_hte_leri.tex"))
}

invisible(mapply(run_leri_models, names(models), models))

# In-party ---------------------------------------------------------------------
party_subgroups <- list(
  reformuk = filter(difdef, inparty == "reformuk"),
  conservative = filter(difdef, inparty == "conservative"),
  labour = filter(difdef, inparty == "labour")
)

party_labels <- c("ReformUK", "Conservative", "Labour")

run_party_models <- function(model_name, formula) {
  
  fits <- lapply(party_subgroups, function(data) lm(formula, data = data))
  
  stargazer(fits[[1]], fits[[2]], fits[[3]],
            title = paste0(toupper(model_name), " for in-party sub-groups"),
            column.labels = party_labels,
            dep.var.caption = "",
            dep.var.labels.include = FALSE,
            align = TRUE,
            omit.stat = c("LL", "ser", "f"),
            no.space = TRUE,
            star.cutoffs = c(0.05, 0.01, 0.001),
            keep = c("type.*", "perception.*", "leri.*", "(Intercept)"),
            notes = "Some variables omitted from table for clarity",
            type = "latex",
            out = paste0("tables/", model_name, "_hte_party.tex"),
            label = paste0("tab:", model_name, "_hte_party.tex"))
}

invisible(mapply(run_party_models, names(party_models), party_models))

# Political news consumption ---------------------------------------------------
news_subgroups <- list(
  low = filter(difdef, pol_news <= 2),
  med = filter(difdef, pol_news == 3),
  high = filter(difdef, pol_news >= 4)
)

news_labels <- c("Low news consumption", "Medium news consumption", "High news consumption")

run_news_models <- function(model_name, formula) {
  
  fits <- lapply(news_subgroups, function(data) lm(formula, data = data))
  
  stargazer(fits[[1]], fits[[2]], fits[[3]],
            title = paste0(toupper(model_name), " for news consumption sub-groups"),
            column.labels = news_labels,
            dep.var.caption = "",
            dep.var.labels.include = FALSE,
            align = TRUE,
            omit.stat = c("LL", "ser", "f"),
            no.space = TRUE,
            star.cutoffs = c(0.05, 0.01, 0.001),
            keep = c("type.*", "perception.*", "leri.*", "(Intercept)"),
            notes = "Some variables omitted from table for clarity",
            type = "latex",
            out = paste0("tables/", model_name, "_hte_news.tex"),
            label = paste0("tab:", model_name, "_hte_news.tex"))
}

invisible(mapply(run_news_models, names(models), models))

# Political interest -----------------------------------------------------------
int_subgroups <- list(
  low = filter(difdef, pol_interest <= 2),
  med = filter(difdef, pol_interest == 3),
  high = filter(difdef, pol_interest >= 4)
)

int_labels <- c("Low political interest", "Medium political interest", "High political interest")

run_interest_models <- function(model_name, formula) {
  
  fits <- lapply(int_subgroups, function(data) lm(formula, data = data))
  
  stargazer(fits[[1]], fits[[2]], fits[[3]],
            title = paste0(toupper(model_name), " for political interest sub-groups"),
            column.labels = int_labels,
            dep.var.caption = "",
            dep.var.labels.include = FALSE,
            align = TRUE,
            omit.stat = c("LL", "ser", "f"),
            no.space = TRUE,
            star.cutoffs = c(0.05, 0.01, 0.001),
            keep = c("type.*", "perception.*", "leri.*", "(Intercept)"),
            notes = "Some variables omitted from table for clarity",
            type = "latex",
            out = paste0("tables/", model_name, "_hte_int.tex"),
            label = paste0("tab:", model_name, "_hte_int.tex"))
}

invisible(mapply(run_interest_models, names(models), models))

# \.\ End of code \.\ #
